Catastrophic cascade of failures in interdependent networks 
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Many systems, ranging from engineering to medical to societal, can only be properly characterized 
by multiple interdependent networks whose normal functioning depends on one another. Failure of 
a fraction of nodes in one network may lead to a failure in another network. This in turn may cause 
further malfunction of additional nodes in the first network and so on. Such a cascade of failures, 
triggered by a failure of a small faction of nodes in only one network, may lead to the complete 
fragmentation of all networks. We introduce a model and an analytical framework for studying 
interdependent networks. We obtain interesting and surprising results that should significantly effect 
the design of robust real-world networks. For two interdependent Erdos-Renyi (ER) networks, we 
find that the critical average degree below which both networks collapse is (k c ) = 2.445, compared 
to (k c ) = 1 for a single ER network. Furthermore, while for a single network a broader degree 
distribution of the network nodes results in higher robustness to random failure, for interdependent 
networks, the broader the distribution is, the more vulnerable the networks become to random 
failure. 



I. INTRODUCTION 

After a decade of intense study on networks, almost all 
work done has concentrated on the limited case of a single 
network which does not interact with other networks [l|, 
U, S, II H 0, 0, B !, M, HO, O • Such situations rarely, if 
ever, occur in nature. Just as in the case of idealized gas, 
when interactions are present as in nature, new physical 
laws appear. 

Analogously, due to technological progress, modern 
systems are becoming more and more coupled together. 
While in the past many networks would provide their 
functionality independently, modern systems depend on 
one another to provide proper functionality. For exam- 
ple, a power network in which the nodes are power sta- 
tions and a communication network in which the nodes 
are computers, are interdependent, since nodes from the 
communication network rely for power supply on the 
power stations, while the power stations depend for their 
control on the proper functioning of the communication 
network. The critical importance of functional depen- 
dence of networks on each other has recently been recog- 
nized 0, G3 • 

In interdependent networks, when nodes in one net- 
work fail, they cause dependent nodes in another network 
to fail. This may happen recursively and can lead to a 
cascade of failures. So a failure of a small faction of nodes 
in only one network, may lead to the complete fragmen- 
tation of all networks. Here, we provide a framework 
for understanding the robustness of interacting networks 
subject to such cascading failures and provide the basic 
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network analytic approach that can underlie future work 
in this area. We present a general model for interde- 
pendent networks that we solve analytically using tools 
from percolation theory and the apparatus of generating 
functions. We present exact analytical solutions for the 
critical fraction of nodes that upon removal will lead to 
a complete fragmentation of all networks. 

Surprisingly, analyzing complex systems as a set of in- 
terdependent networks may destabilize the most basic 
assumptions that network theory has relied on for single 
networks. While for a single network a broader degree 
distribution of the network nodes results in the network 
being more robust to random failure, for interdependent 
networks, the broader the distribution is, the more vul- 
nerable the networks become to random failure. The 
implications are dramatic - the current methods applied 
to the design of robust networks need to be modified to 
include the findings about interdependent networks. 



II. THE MODEL 

Consider two networks A and B and assume that the 
functioning of a node in network A depends on the ability 
of one or more nodes in network B to supply a critical 
resource to the node in network A. Similarly, a node in 
network B depends on a set of nodes in network A. The 
networks can be connected in different ways; in the most 
general configuration one could specify the distributions 
of connections between the nodes from both networks. 

The networks can have the same, or different, topolo- 
gies. The model can easily be extended to an arbitrary 
number of interacting networks each with its own specific 
topology and dependence on the other networks. For ex- 
ample, an interesting dependence for three interacting 
networks could be a circular dependency in which the 
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nodes in network B depend on network A for a resource, 
the nodes of network C depend on the nodes of network 
B for a resource and the nodes of network A depend on 
network C for resources. 

Our key insight from percolation theory is that for each 
of the networks to remain functional after nodes have 
failed, the network must include a spanning cluster of 
functional nodes. Nodes that are not part of the span- 
ning cluster will become nonfunctional and will cause the 
nodes from other networks that are connected to them 
to also become nonfunctional. 

For simplicity and without loss of generality, we will 
assume a system of two networks, A and B, both with 
N nodes. Within each network, we assume that nodes 
are randomly and independently connected with degree 
distributions Pa (A) and Ps(fc). We also assume, for sim- 
plicity, that each node in network A is connected to, and 
dependent on, one node in network B and vice-versa. 
Next we will remove a fraction of nodes 1 — p from net- 
work A and all the edges connected to these nodes, so 
that only a fraction p remain functional. Simultaneously, 
the corresponding nodes (and their edges) in network B 
are removed since they are dependent on the nodes in 
network A. 

As edges are removed, the networks break up into con- 
nected components (clusters). The clusters in network 
A and the clusters in network B will be different be- 
cause the networks are each connected differently We 
define a mutually-connected cluster as the set of nodes 
in network A which belong to a cluster in network A 
and also have their corresponding nodes in network B 
belong to a single cluster in network B (or vice- versa). 
We assume that clusters of nodes that are disconnected 
from the network core (giant component/spanning clus- 
ter) become non-functional and are removed. Thus, the 
mutually-connected giant component will be of special 
interest since it is the only functional part of the system. 

The questions that we will ask: What is the critical 
p = p c below which all the mutual clusters constitute only 
an infinitesimal fraction of the network, i.e., no mutual 
giant component exist? What is the probability Paoip) 
for a node to belong to the mutual giant component as 
function of pi To solve this model we will introduce a 
recursive process which we will identify with a physically 
meaningful cascade of failures. 

To solve this model we will first define the ai-clusters 
of network A after only a fraction of nodes p remain. 
Next we will treat each of these ai-clusters as separate 
subsets of a network B, i.e. all the B-edges connecting 
different ai-clusters will be removed. We will define this 
state of the networks as the first stage in the cascade of 
failures. Accordingly, each of the ai-clusters may split 
into several 62-clusters. Some of the ai-clusters will not 
split and will coincide with b 2 clusters. Obviously such 
clusters are mutually connected. Finally we remove from 
the network all A-edges connecting different 62 clusters. 
We will define this state of the networks as the second 
stage in the cascade of failures. Analogously, in the third 



stage we will determine all the a3-clusters and in the 
fourth stage we will determine all the 64-clusters, and 
will continue this process until no further splitting and 
edge-removal will occur. 

Note that in this process the majority of new mu- 
tual clusters identified after each stage of failures will 
be isolated nodes, few of them will be of size 2 and very 
rarely we will have larger mutual clusters. Indeed if we 
have two nodes that are connected by an A-edge, the 
probability that they will be connected by a B edge is 
1 - E k PB(k)(l - l/N) k » 1 - £ fc P B (fc)(l - k/N) = 
T, k p B(k)k/N =< k > B /N — ► for N — ► 00. The 
probability that three nodes connected by A-edges are 
also connected by B edges scale as 1 /N 2 and so on. 



III. ANALYTICAL SOLUTION 

Now we will solve the problem analytically using the 
apparatus of generating functions. As in Refs. [l5l [l6j 
we will introduce generating functions of the degree dis- 
tributions 



and 



G A0 (x) = J2 p A(k)x k 



G so (x) =]TP s (fc)s 



(1) 



(2) 



Analogously we will introduce generating functions of the 
underlining branching processes: 



and 



G A1 (x)=G' A0 (x)/G' A0 (l) 



G m (x) = G B0 (x)/G' B0 (l). 



(3) 



(4) 



Random removal of fraction 1 — p of nodes will change 
the degree distribution of the remaining nodes [15j], so 
that the new generating functions become 



and 



Gao(x,p) = G A o(l -p(l - x)), 



Gbo(x,p) = G B o(l - p(l - x)), 



G A i(x,p) = G A i(l -p(l - x)), 



G B1 (x,p)=G B1 (l-p(l-x)). 



(5) 



(6) 



(7) 



(8) 



Let us denote the subset of nodes after the random re- 
moval of 1 — p nodes a,s Aq — Bq C A = B . If the number 
of nodes in the entire network is N, the number of nodes 
in Ao — Bq is Nq = pN. The fraction of nodes that 
belong to the giant component of the network Aq is 



Pa(p) = 1 - Gao(/a,p), 



(9) 
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where f A satisfies a transcendental equation 

fA(p)=GAl(fA,p). (10) 

Equation pop can be simplified by substitution z A — 
l-p(l-/ x ) 

l-l/p + z A /;P =6.41(2.4). (11) 
Then Eq. © becomes 

p A (p) = l-G A0 (z A ), (12) 

Analogous equations characterize the giant component 
of network Bq. After the initial attack which removes 
(1 — p)-fraction of nodes from both networks, the first- 
stage failure is caused by the fragmentation of the subset 
A . The giant component Ai of A will constitute P A (p)- 
fraction of Aq. Thus the number of nodes in A\ is N\ = 
NqP a { P ) =pP A (p)N =p x N. 

After the first-stage failure the fraction of functioning 
nodes is p\ = pP A {p) (subset Ax). Because the nodes of 
the networks B and A coincide, the same fraction of nodes 
remains functioning in network B. Because the topology 
of network B is independent the topology of network A, 
these functioning nodes are totally random with respect 
to connections in network B. Thus we can again apply the 
apparatus of generating functions and find the fraction 
Pb(pi) of the giant component B>2 of network B with 
respect to the subset A\. The number of nodes in the 
giant component B 2 C Ai is N 2 = p 2 N = P B {pi)Nx = 
P B (pi)piN = pP A {p)P B \pi)N . Thus, the fraction of 
functioning nodes after the second stage failure is p 2 — 
pPa(p)Pb(pi) (subset B 2 ). 

Now we will analyze what happens during the third- 
stage failure which is caused by further fragmentation 
of the giant component A\ by removal N\ ~ N 2 — 
(1 — Pb(pi)Ni, nodes which do not belong to B 2 . The 
removal of these nodes form A\ is equivalent to the re- 
moval of the same fraction of nodes from Aq (because 
all the nodes that were removed at the stage of the ini- 
tial attack do not belong to B 2 , A\, and Aq. The total 
number of nodes that must be removed from network 
A is (1 — Pb{pi))Nq nodes from Aq plus the number of 
the initially attacked nodes (1 — p)N. Thus, the total 
number of nodes that must be removed from network 
A is (1 — pPb{pi))N. Hence the third-stage failure is 
equivalent to a random attack in which p is replaced by 
p 2 = pPb{pi). Accordingly the number of nodes in the 
giant component A 3 C B 2 is N 3 = p s N — p' 2 P A (p' 2 ). 

Following this approach we can construct the se- 
quence of giant components in the cascade of failures: 
A 2m +i c B 2rn c A 2rn -i . . . c A 3 c B 2 c Ai c 
Bq = Aq C A = B. The number of modes in 
each giant component of this sequence is N > pN = 
NpQ > Npi > ... Np 2m+ i . . ., where the numbers p n 
can be obtained by recursive relations: pq = p' Q = p, 
Pi = Pi = PqPa(p'q), p'i = pPb(p[), P2 = p[Pb(p'i), 
Ps = pPa(p' 2 ), P3 = P 2 Pa(p' 2 ) ■ • -p' 2m = pPB{p' 2m -i), 



P2m = P2m-l P B{P2m-l), P2 m +1 = P P A{p' 2m ), P2m+1 = 
P2m P A{p' 2m )- 

Now we will determine the size of the mutual giant 
component. The fraction of nodes in the mutual giant 
component, Poo is the limit of the sequence p n for n — > 00. 
This limit must satisfy the equations p 2rn +i — P2m — 
p 2m -i since the cluster is not further fragmented. Using 
relations between p n and p' n _i, and denoting p' 2m ^i = x 
and p' 2m — y we arrive to a system of two symmetric 
equations with two unknowns: 

(x=pp A (y) , 13 
\V=PPb{x). 

This system of equations has one trivial solution x — 0, 
y = for any p, corresponding to the zero size of the 
giant mutual component. If p is large enough there exists 
a different solution which gives the nonzero size of the 
mutual giant component. We can easily exclude y from 
these equations and obtain a single equation 

x = pp A (pp B (x)) (14) 

This equation can be solved graphically (Fig[2]) as the 
intersection of a straight line y = x and a curve y = 
PPa(pPb(x)) which both intersect at the origin. When 
p is small enough the curve increases very slowly and 
does not intersect with the straight line. The criti- 
cal case when the nontrivial solution emerges, corre- 
sponds to the case when the line touches the curve at 
a single point x and in this point we have a condition 
1 = p 2 p' A (jppB{x))p' A {x), which together with equation 
x = pp a (ppb(x)) gives the solution for the critical p and 
the critical size of the mutual giant component. 



IV. ER NETWORKS 

In case of ER networks, whose degrees are Poisson- 
distributed (l7l [l8j , the problem can be solved explicitly. 
Suppose that the average degree of the network A is a and 
the average degree of the network B is b. Then, G A \ (x) = 
G A q = exp(a(x — 1)) and Gbi = Gbo = exp(6(a; — 1). 
Accordingly system (fT3"|) becomes 

l X=P \]- f f A } (15) 

\y = p[l-fB\, 

where 

f f A = exp[ay(f A - 1)] , , 

\ f B = exp[bx(f B - 1)]. [W) 

Excluding x and y, we get a system with respect to f A 
and Jb ■ 

( f A = e -ap(/ A -l)(/B-D 
\ f B = e -&P(/A-l)(/B-l). 
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Introducing a new variable r = fV a — f^ b , we reduce 
system (jTTJ) to a single equation 



r = e -p(r a ~l)(r b -l) 



(18) 



which can be solved graphically for any p. The critical 
case corresponds to the tangential condition 



The right-hand side of this equation behaves as x^ , where 
/i = l/[(3 — A^)(3— Xb)] > 1. Thus the r.h.s curve always 
goes below y — x for x — > 0, so for sufhciently small p 
we do not have a non-trivial solution, which means the 
absence of the mutual giant component. Thus we have a 
percolation transition at some p = p c > 0. 



1 = JL e -P(r a -l)<r -1) = p [ar a + br b -{a + b)r a+b }, (19) 
dr 

from where the critical value of r = r c satisfies transcen- 
dental equation 



— 



(20) 



and the critical value of p = p c can be found from 
Eq. IH]). 



1 



Pc 



ar^. + br b c — (a + b)r, 



a+b ' 



(21) 



The values of p c and Px, for different a and b > a are 
presented in Fig. [3] as function of a/6. In case a = b, 
Ia = /fl = /j and / c satisfy equation 



d-/c) 2 

/ c = e 2^+ 2 / c 



(22) 



which gives a solution / c = 0.28467, p c = 2.4554/a, and 
the critical fraction of nodes in the mutual giant compo- 
nent Poo = p c (l — fc) 2 = 1.2564/a. Numerical simula- 
tions of the ER networks are in excellent agreement with 
the theory (Fig. ij. 



V. SCALE-FREE NETWORKS 

For regular percolation in a scale-free network with a 
power law degree distribution Pa (A;) ~ kr XA , it is known 
that p c — > 0, as N — > oo for Xa < 3. Surprisingly, for mu- 
tual percolation this is not the case and p c remains finite 
for Xa > 2. To see this, we can find analytical approx- 
imation for Pa(p)- First, we begin by solving Eq. pTjl . 
According to Tauberian theorems, for A^ < 3, Gai{x) 
has a singularity at x = 1 of the sort 1 — «u(l — x) Xa ~ 2 . 
Therefore it has a diverging derivative which has a phys- 
ical meaning of the branching factor Ha- To solve Eq. 
(|lip . we must find the intersection of the straight line 
y = 1 — 1/p + z/p and the curve y = Gai(z). The 
straight line passes through the point y = 1, z = 1 with 
the derivative 1 /p. Thus there is always a trivial solution 
z = 1, which corresponds to the absence of percolation. 
If, G[ 41 (l) = is finite, we do not have another solution 
for p < 1/kA (a classical result for regular percolation), 
but for Xa < 3, we always have a non trivial solution 
z = 1-{pka) 1/( - 3 ~ Xa} - Since G' A0 (1) = (fc), which is finite 
for X A > 2, Eq. (TfJ) yields P A (p) = (pK A ) 1/(3 ~ Aj4) ( fc ^> ■ 
Finally Eq.(fT4|) becomes 



z = p(fcA> [p^bX^) 1 ^ 3 -^! 1/(3 . (23) 



VI. ROBUSTNESS OF INTERDEPENDENT 
NETWORKS 

For interdependent networks we find the surprising be- 
havior that networks with a broad degree distribution (of 
the network nodes) are more vulnerable to random at- 
tack compared to networks with a narrow distribution. 
To understand this result we note the following: 1) All 
interdependent networks are randomly connected, high 
degree nodes from one network might be connected to 
low degree nodes from the other networks. 2) At each 
step when nodes (and their links) are disconnected from 
one network their corresponding nodes (and their edges) 
from the other network are also removed. 

Therefore, the hubs that play such a dominant role 
in the robustness of single networks become vulnerable 
when a cascade of failures occurs in interdependent net- 
works. Moreover, for a network with a fixed average de- 
gree, a broader distribution means more nodes with low 
degree to balance the high degree nodes. Since the low 
degree nodes are more easily disconnected the advantage 
of a broad distribution on single networks becomes a dis- 
advantage for several interdependent networks. 

In Fig. [7] we compare simulation results for several SF 
networks with different A values, an ER network and a 
Random Regular (RR) network, all with an average de- 
gree (fc) = 4. The simulation results are in full agreement 
with our analytical results and it can clearly be seen that 
for a broader distribution p c is indeed higher. 



VII. FINITE SIZE EFFECTS 

Our considerations are rigorous for N — * oo. For a 
finite network, the relative fluctuations of all fractions 
decrease as lVN so, for the finite network, there is a 
range of values of p for which the mutual giant compo- 
nent exists with probability Pqo(p) (Fig. [5]). Its deriva- 



tive diverges as N — > oo as dP oa /dp\ p=Pa ~ \N, and 
for TV — * oo, Poo(p) becomes a step function Poa{p) = 
for p < p c and P OQ {p) — 1 for p > p c . The square-root 
scaling with N of the width of the interval p for which 
we can have a complete fragmentation for some realiza- 
tions of networks and a giant component for the other 
realizations of the networks can be justified by the fol- 
lowing arguments. The actual fraction of the remaining 
nodes p a in a finite network of size iV will be normally 
distributed around given p with the standard deviation 
inverse proportional to N. Thus Poo(p) is equal to the 
probability that p a > p c , which is equal to the integral of 



5 



the normal probability density with zero mean and the 
same standard deviation from p c — p to infinity. There- 
fore the derivative dP^/dp has a Gaussian shape with 
standard deviation proportional to 1/y/N. 

The average number of stages (n) in a cascade of fail- 
ures for p > p c diverges proportionally to InN/^/p — p c . 
This follows from the properties of the iterative process 
Eq. This can be seen from the fact that near 

p = Pc Eq. (fTi)) has two roots produced by the inter- 
section of the curvy line which can be approximated by 
a parabola y = a(p)x 2 + b(p)x + c(p) and a straight line 
y = x (Fig. [2]). This is equivalent to solving a quadratic 
equation a{p)x 2 + (b(p) — l)x + c(p) — 0. The value 
p = p c is given by the discriminant of this equation 
equal to zero: d(p c ) = {b{p c ) - l) 2 - Aa(p c )b(p c ) = 0. 
In the general case, all three parameters, a(p), b(j>), and 
c(p), have non-zero derivatives at p — p c . Therefore, 
in the general case d{p) has also a non-zero derivative 
at p — Pc, and hence the difference between the roots 
scales as s/p — p c . Thus, the derivative of the curve at 
the largest root, which corresponds to the limit of the 
iterative process scales as /' = 1 — U\Jp — p c , where a is 
some positive constant. For Eq. (fl~3|) the iterations con- 
verge to the root as f' n = exp(— cty/p — p c n). In a real 
network, they will stop when the difference between two 
successive iterations will be smaller than one node, which 
yields a condition exp(— cty/p — p c n) ~ 1/N. Hence in- 
deed (n) ~ In N/ ^Jp — Pc- 



For p < p c the solution does not exist and the curve 
misses the line with the distance proportional to the neg- 
ative descriminant. As the curve comes close to the line 
the steps are proportional to (x— x c ) 2 +d, where d ~ p c —p 
is the minimal distance between the curve and the line. 
The number of such steps per dx is dx/((x — x c ) 2 + d). 
The total number of steps are thus the integral of this 
quantity between x = p and x = 0, which in the limit 
d — > Q gives (n) = tt / ' \fd ~ l/y/p c — p. 

Exactly at the critical point p = p c the straight line 
touches the curve at a single point and the sequence of 
iterations converges as x n +\ — x c = x n — x c — a(x n — x c ) 2 . 
These iterations converge to x c as 1 jn which can be seen 
by plugging into this equation x n — x c — C /n^ + o(n~' 3 ) 
where C and (3 are some unknown constants. Expanding 
(ri + I)-! 3 in Tailor series and equating coefficients for 
equal powers, one can see that (3=1. However, in real 
network, due to Gaussian spread in p a , we are never at 
criticality, and the typical p a — p c ~ ly/~N. Therefore the 
distributions of the number of stages in the cascade has 
an exponential tail exp[— an^/p a — p c ], in which {p a —Pc) 
must be replaced by its typical value 1/^/N. Therefore, 
the distribution of P(n) must have an exponential tail 
P(n) ~ exp[— a'n/N 1 / 4 ], where a' is some positive con- 
stant. Thus at criticality, we expect that (n) ~ N 1 / 4 as 
supported by our simulations (Fig. [S]). 
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FIG. 1: (Color online) Demonstration of two interdependent networks. Nodes in network B (communications network) are 
dependent on nodes in network A (power grid) for power; nodes in network A are dependent on network b for control information. 
General case is represented in which there is not a one-one correspondence between nodes in networks. 




FIG. 2: (Color online) Iterative process described by Eq. (|14[) for the case of the scale-free distribution Pa(/d) = Ps(fc) = 
(2/fc) 2 — (2/(fe + l)) 2 for k — 2,3,.... For k — > oo, this distribution scales as fc _A , where A = Aa = As = 3. Three curves 
corresponding to p — 0.70 < p c (black), p — 0.752 « p c (red) and p — 0.80 > p c (green). One can see that for p > p c , the 
iterations (red arrows) starting from p' — p, converge to the largest of the two roots of Eg . (I14fl . For p < p c , the iterations 
converge to 0. 




FIG. 3: ER networks: critical fraction p c and the fraction of nodes in the mutual giant component at criticality as function 
of the ratio a/b, where a and b are the average degrees of networks A and B respectively. 
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On 



— N=128000 
G-© theory 




FIG. 4: (Color online) Comparison of the fraction of the giant components after n stages of the cascade failures for several 
random realizations of ER networks with a = b, N — 128000 and ap — 2.45 < ap c = 2.4554 and theoretical prediction of 
Eq. (|13[l , One can see that for the initial stages the agreement is perfect, however at larger n the deviations due to random 
fluctuations of the order of 1/yN in the actual fraction of the remaining nodes p a start to increase. The theoretical prediction 
after a region of the plateau around the the critical value, drops to zero, corresponding to the complete fragmentation of the 
network. The random realizations separate into two classes: one that converge to a mutual giant component and the other 
that results in a complete fragmentation. 
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FIG. 5: Numerical simulations of ER networks with a = b and finite number of nodes, N. The probability of existence of the 
mutual giant component Poo, is shown as function of p for different N. One can see that as N — > oo the curves converge to a 
step function. The theoretical prediction of p c is shown by the arrow. 




FIG. 6: Scaled distribution of the number of stages in the cascade failures for ER graphs with a = b at criticality (pa 
for different values of N. 
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FIG. 7: (Color online) Simulation results for as a function of p for for SF networks with A = 3, 2.7, 2.3, an ER network and 
a Random Regular (RR) network, all with an average degree (k) = 4. The simulation results are with full agreement with our 
analytical results and it can clearly be seen that for a broader distribution p c is higher. 



